Land Cover Classification¶
Raster maps are geo-reference images, in which infomration is classified in the pixel values. In simple terms, the difference between regular images and geo-referenced raster images is that the latter has each pixel associated to a geographic (or projected) coordinate. So instead of just using rows and columns like we would do in regular images, we can access data in raster images using latitude and longitude.
The following example uses a land cover map for 2019 for Barton county, KS. This county is located in central Kansas and has a bit of everything: towns, lakes, roads, cropland, and grassland. The goal of this exercise is to learn how read raster data, access the individual values and coordinates, apply some operations to the data to extract features, and use image analysis to extract features of potential interest.
import numpy as np
import holoviews as hv
from holoviews import opts
hv.extension('bokeh')
import xarray as xr
# Load maps
B2010 = xr.open_rasterio('../datasets/land_cover/barton_county_ks_2010.tif')
B2019 = xr.open_rasterio('../datasets/land_cover/barton_county_ks_2019.tif')
# Show map
options = opts.Image(cmap='Spectral',aspect='equal',colorbar=True)
barton_2010_map = hv.Image( (B2010.x,B2010.y,B2010.values[0]),label='Barton County, KS 2010').opts(options)
barton_2019_map = hv.Image( (B2019.x,B2019.y,B2019.values[0]),label='Barton County, KS 2019' ).opts(options)
(barton_2019_map).opts(frame_width=300) # Use this to dsiplay both figure: (barton_2019_map + barton_2019_map).opts(frame_width=800).cols(2)
Map codes¶
The image has codes for each land cover encoded in uint8, which means unsigned integers of 8 bits. The values range from 0 to 255, meaning that this way of encoding information can only store 256 possible land covers, which is plenty for most applications. For example a pixel with a value of 111 representes areas occupied by water in 2019 and pixels with a value of 24 represents fields planted with winter wheat in 2019. The total of 256 land covers comes from the fact that a two-bit system would allow us to represent only two land covers using either 0 or 1. We actually use this format quite often when we plot boolean arrays.
To expand the number of possibilities we need more slots with 0s and 1s. So, since \(2^8=256\), using eight digits that can be either 0s or 1s allows for a total of 256 combinations. So, the value for grassland 176 would be represented in binary as 10110000 You can try this on your own by running the following Python code: bin(176). The first two values of the printed sequence (something like 0b) simply mean to tell you that the numbers are represented using the binary system.
Land Cover Codes¶
Land cover |
Code |
Comments |
|---|---|---|
Background |
0 |
|
Corn |
1 |
|
Cotton |
2 |
|
Rice |
3 |
|
Sorghum |
4 |
|
Soybeans |
5 |
|
Wheat |
24 |
Winter wheat and Spring wheat |
Open water |
111 |
Lakes, ponds, rivers, streams and creeks |
Urban areas |
121, 122, 123, 124 |
Each code is for very low to high development area |
Barren land |
131 |
|
Forests |
141, 142, 143 |
Refers to deciduous, evergreen and wetland forests |
Shrublands |
152 |
|
Grassland/Pastures |
176 |
|
Wetlands |
190 |
# Print unique map codes
np.unique(B2019.values)
array([ 0, 1, 2, 4, 5, 6, 21, 24, 26, 27, 28, 29, 31,
36, 37, 44, 53, 58, 59, 61, 111, 121, 122, 123, 124, 131,
141, 142, 143, 152, 176, 190, 195, 205, 225, 228, 236, 238],
dtype=uint8)
# Classify some land covers
water_2019 = B2019 == 111
urban_2019 = (B2019>=121) & (B2019<=124)
wheat_2019 = B2019 == 24
grassland_2019 = B2019 == 176
# Plot land covers
water_map_2019 = hv.Image( (water_2019.x,water_2019.y,water_2019.values[0]) ).opts(cmap='Blues',title='Water bodies')
urban_map_2019 = hv.Image( (urban_2019.x,urban_2019.y,urban_2019.values[0]) ).opts(cmap='Purples',title='Urban area')
wheat_map_2019 = hv.Image( (wheat_2019.x,wheat_2019.y,wheat_2019.values[0]) ).opts(cmap='Reds',title='Wheat')
grassland_map_2019 = hv.Image( (grassland_2019.x,grassland_2019.y,grassland_2019.values[0]) ).opts(cmap='Greens',title='Grassland/Pasture')
(water_map_2019 + urban_map_2019 + wheat_map_2019 + grassland_map_2019).cols(2)